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A  nonlinear  progressive  wave  equation  (NPE)  describes  the  evolution  of  a  low  Mach  number 

shock  wave.  The  NPE  is  the  nonlinear  time  domain  counterpart  of  the  frequency  domain 
linear  parabolic  wave  equation  (PE)  for  small  angle  propagation.  The  NPE  in  spherical 
symmetry  admits  a  similarity  solution  that  specifies  both  the  shape  of  the  pulse  and  the  shock 
strength  as  a  function  of  range.  For  finite  amplitude  spherical  waves  whether  self-similar  or 
not,  the  theory  predicts  constancy  of  an  impulse  integral  corresponding  to  that  of  linear  theory. 
For  the  self-similar  waves,  theory  and  available  data  are  in  qualitative  agreement  in  the 
following  areas:  ( 1 )  The  shock  strength  decreases  with  range  as  an  approximate  power  law; 

(2)  the  temporal  behavior  of  the  solution  at  fixed  range  is  a  shock  discontinuity  followed  by  a 
roughly  exponential  decay;  and  ( 3 )  the  effective  relaxation  time  behind  the  shock  is  in  ^ 

reasonable  agreement  with  data  for  slightly  more  than  the  first  decade  in  range.  \Ji-f '  <  ' 
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INTRODUCTION 

A  number  of  sources  give  rise  to  waves  that  spend  much 
or  all  of  their  lifetimes  as  “weak”  shocks,  propagating  only 
slightly  faster  than  Mach  1.  These  sources  include  super¬ 
sonic  aircraft,'  underwater  explosions,2-3  and  a  number  of 
laboratory  experiments  and  procedures.  A  shock  is  consid¬ 
ered  mathematically  weak  when  the  density  of  the  medium 
changes  by  a  small  amount  (say,  20%  or  less)  across  the 
shock  front.  Underwater  shocks  are  almost  always  weak, 
even  though  overpressures  may  be  as  high  as  thousands  of 
atmospheres.  If  overpressure  is  used  to  measure  the  strength 
of  the  shock,  it  is  important  to  keep  in  mind  that,  for  water, 
the  relevant  scale  for  comparison  is  not  ambient  pressure  as 
in  the  case  of  a  gas,  but  bulk  modulus  -  p<f\ ~22  X  103  atm. 
Here,  p0  and  c0  are  the  density  and  linear  sound  speed  of 
water. 

Experimental  data  from  underwater  detonation  of  high- 
explosive  charges  reveal  an  empirical  power  law  dependence 
of  the  shock  strength  (i.e.,  overpressure)  versus  range  from 
the  source.4  Pressure  relaxation  times  tp  =  \p  (dp/dt)~'\ 
behind  the  shock  front  also  fall  with  somewhat  more  scatter 
about  an  approximate  power  law  in  range.  Specifically,  peak 
overpressures  and  relaxation  times  for  underwater  explo¬ 
sions  vary  with  a  scaled  range  coordinate  according  to 

pm.%  =0.024Poc2(»Fl/3/r)"3, 

0.4 5 r/ W' °  5400  m  kg~'/3, 
rp  =92.5xl0-6IF,/3(r/IFl/3)022,  (I) 

0.4Sr/IF,/3S40mkg-''3, 

where  r  is  the  distance  from  the  explosion  in  meters,  W  is  the 
TNT  equivalent  mass  of  the  explosive  charge  in  kilograms, 
and  rp  is  in  seconds.  Overpressures  from  a  number  of  shots 
involving  different  explosive  charges  fall  within  about 
±  20%  of  the  power  law  over  approximately  3  decades  in 
the  scaled  range.2  Relaxation  time  data  are  subject  to  more 


scatter  than  the  overpressure  data.  The  scatter  about  the 
power  law  ( 1 )  for  relaxation  time  is  approximately  ±  30%. 

The  scaling  of  range  by  W 1 13  in  ( I )  is  easily  explained  by 
dimensionality  considerations.  When  a  point  source  dumps 
hydromechanical  energy  E  impulsively  and  symmetrically 
into  an  infinite  homogeneous  medium  whose  pressure  is  neg¬ 
ligible,  the  only  macroscopic  scale  size  inherent  to  the  con¬ 
figuration  is  L0  =  ( E /ptf  o ) 1  /3,  referred  to  by  Sedov5  as  the 
dynamic  range.  Here,  L0  is  an  order-of-magnitude  estimate 
for  the  shock  radius  at  a  time  when  the  root-mean-squared 
fluid  particle  speed  falls  below  e0.  The  energy  E  is  propor¬ 
tional  to  the  mass  W  of  the  explosive,  so  that  scaling  to  L0 
implies  scaling  to  W'n.  The  scaling  of  rp  by  fVtn  reflects  the 
inherent  time  scale  Lo/c0. 

The  approximate  power  law  behavior  of  the  overpres¬ 
sure  and  relaxation  time  with  range  is  not  so  obvious.  The 
simplicity  of  the  range  dependence  of  data  summarized  by 
( 1 )  does,  however,  suggest  some  sort  of  similarity  behavior 
in  the  evolution  of  the  shock  wave.  Similarity  solutions  have 
long  been  known  to  be  relevant  to  strong  shocks,5  6  but  as¬ 
sumptions  regarding  the  dominance  of  the  shock  over  the 
ambient  medium  fail  as  the  shock  grows  weak. 

Rogers2  presented  a  semiempirical  theory  whose  nu¬ 
merical  solution  reproduced  approximately  the  r  ~  1 1 3  power 
law  for  overpressure  versus  range.  The  range  behavior  of  the 
relaxation  time  from  this  theory  was,  however,  less  satisfy¬ 
ing.  The  predicted  relaxation  time  increased  too  rapidly  with 
range  close  to  the  source  and  then  too  slowly  far  from  the 
source.  The  empiricism  in  this  method  lay  in  imposing  the 
assumption  (with  some  support  from  available  data)  that, 
just  after  the  passage  of  the  shock  past  a  fixed  range  point, 
the  temporal  behavior  of  the  overpressure  was  an  exponen¬ 
tial  falloff.  Actually,  this  behavior  can  hold  only  for  a  limited 
time,  for  the  density  profile  must  change  sign  at  some  point 
in  order  to  conserve  mass.  This  will  be  discussed  further 
below.  The  solution  method  of  Ref.  2  employed  a  one-di- 


mensioi^l  njetbod  of  characteristics  immediately  behind  the 
shock,  followecf*by  an  approximate  adaptation  to  spherical 
Symmetry.  $>~\1 

In  contrast,  the  metbrfd  used  here  begins  with  a  reduced 
nonlinear  wave  equatkmaccurate  to  second  order  in  the  den- 
•  s»ty  perturbation.  A  fq^jnafsimilarity  solution  is  extracted  in 
'  which  time  and  rangl?4tependence  enter  through  a  single 
variable  Our  results  show  a  more  rapid  dropoff  of 
overpressures  near  the  source  than  Rogers’  predictions.  The 
relaxation  time  prediction  from  our  theory  is  in  better  agree¬ 
ment  with  data  at  close  ranges,  but  also  flattens  off  far  from 
the  source. 

I.  THE  NONUNEAR  PROGRESSIVE  WAVE  EQUATION 
(NPE) 

The  propagation  of  three-dimensional  finite  amplitude 
acoustic  pulses  subject  to  refraction,  diffraction,  and  nonlin¬ 
ear  steepening  is  described  by  a  nonlinear  progressive  wave 
equation  (NPE)  derived  in  a  related  article.7  The  NPE  has 
been  shown1  to  be  the  nonlinear  time  domain  equivalent  of 
the  frequency  domain  parabolic  wave  equation  (PE).9  The 
NPE  for  a  spherical  nonlinear  wave  in  a  homogeneous  medi¬ 
um  is 

where 

X  ~p\/po -  (3) 

Here,  px  is  the  density  perturbation,  and 

+  -,+lrl£l  <4> 

2cq  op  o  o  log  p  I  o 

InEq.  (4)0  is  a  thermodynamic  quantity  referred  to  as  the 
medium’s  coefficient  of  nonlinearity.  For  a  perfect  gas  with 
specific  heat  ratio  y,  0  =  ( 1  +  y)/2.  The  value  of  (3  for  wa¬ 
ter  is  approximately  3.5.  The  derivatives  in  (4)  are  taken 
from  an  adiabatic  equation  of  state  for  the  medium.  This  is 
permissable  in  our  model,  since  for  a  weak  shock,  heating 
(entropy  increase  as  opposed  to  adiabatic  compression)  is 
proportional  to  the  cube  of  the  density  jump. 10  An  estimate 
of  the  smallness  of  shock  heating  can  be  obtained  from  Table 
2. 1  of  Cole3  for  a  water  shock  of  overpressure  Ap  =  5000 
atm  (an  approximate  upper  limit  to  the  region  of  applicabili¬ 
ty  of  our  model).  Immediately  behind  the  shock,  the  total 
temperature  increase  (adiabatic  compression  plus  entropy 
generated  at  the  shock)  A  7" is  approximately  S  *C.  Values  for 
T  AS  reveal  that,  upon  return  to  atmospheric  pressure,  the 
water  is  heated  by  only  about  1.3  *C.  We  maintain  that  this 
reflects  comparison  of  a  cubic  to  a  quadratic  effect  when  Ap/ 
p<f o  =0.25.  (The  fact  that  the  constant  of  proportionality  in 
this  comparison  is  near  unity  is  fortuitous. ) 

Equation  (2)  describes  the  evolution  of  a  finite  ampli¬ 
tude  compressional  wave  in  the  farfield  of  a  spherical  source. 
Inherent  in  the  derivation  of  the  nonlinear  term  in  (2)  is  the 
assumption  that  it  applies  to  a  wave  whose  gradient  scale  size 
at  a  point  is  much  less  than  the  radius  r  from  the  origin.  In 
the  absence  of  the  nonlinear  term,  ( 2 )  describes  exactly  out¬ 
going  spherical  wave  solutions 

R=Ar-Cot)/r  (5) 
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of  the  linear  wave  equation  in  spherical  symmetry: 

^~  =  cl72R.  (6) 

In  this  solution,  mass  conservation  relies  on  a  consistency 
condition  referred  to  as  the  conservation  of  impulse: 

J /U)<&=J  rR  dr  —  0.  (7) 

The  integration  limits  in  (7)  extend  over  the  entire  range 
where/is  nonzero.  The  second  form  of  this  condition  will  be 
shown  also  to  apply  to  nonlinear  spherical  waves.  Since  Eq. 
(2)  retains  only  linear  and  quadratic  terms,  the  assumption 
of  constant  (3  allows  (3  to  be  absorbed  into  a  scaled  density 
perturbation 


r  =  #>,/p0- 

The  resulting  equation  for  s  is  then 

(8) 

w 

& 

+ 

T" 

i 

it 

•SI* 

(9) 

where 

<!>(S)  =  s  +  sV2. 

An  immediate  result  of  (9)  is  that  the  location  of  a  contin¬ 
uous  zero  crossing  (s  =  0)  propagates  at  a  speed  c„. 


It.  SIMILARITY  SOLUTION 


A  similarity  solution  to  (9)  may  be  found  by  combining 
the  independent  variables  r  and  t  into  a  single  independent 
variable  Cgt  /r.  A  solution  is  sought  in  the  form 


s  =  r(V).  q-cj/r. 

(10) 

Substituting  ( 10)  into  (9)  leads  to 

no 

The  nonlinear  partial  differential  equation  (9)  is  thus  con¬ 
verted  to  an  ordinary  differential  equation  for  a  restricted 
class  of  solutions  of  the  form  ( 10) .  Equation  (11)  may  be 
integrated  for  general  by  interchange  of  dependent  and 
independent  variables  to  yield 

vm 1121 

For  ^(s)  as  given  in  (9),  (12)  becomes 

V-  1  +s  +  s(l  +  s/2)log[|eu|/(s  +  2)],  (13) 

where  a  is  a  constant  of  integration.  Equation  (13)  implicit¬ 
ly  defines  5  as  a  function  of  rj,  i.e.,  the  self-similar  overdensity 
profile  as  a  function  of  range  and  time.  Equation  (13)  gives 
i)-»  1  as  r-*0,  leading  to  a  zero  crossing  at  r  =  c„r  for  all  a. 

III.  SHOCK  PROPAGATION 

We  seek  solutions  describing  shock  jumps  with  an  un¬ 
disturbed  medium  ahead  of  the  discontinuity  and  the  contin¬ 
uous  profile  (13)  behind  it  (or  between  discontinuities  in  the 
case  of  an  N  wave  shock).  We  must  show  that  the  jump 
conditions  at  the  shock  can  be  satisfied  by  the  similarity  vari¬ 
ables  in  a  consistent  way.  The  form  of  Eq.  (9)  is  such  that 
results  of  shock  physics  can  be  applied  directly.  Assuming  a 
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discontinuity  between  values  to  the  left  and  right  of  the 
shock  (denoted  by  subscripts  /  and  r),  the  requirement  of 
mass  conservation  across  the  shock  specifies  the  speed  at 
which  the  shock  front  propagates  ( i.e.,  the  Rankine-Hugon- 
iot  relation): 


dt  s,  —  s,  \  2) 


The  second  expression  in  ( 14)  assumes  that  the  shock  is  a 
jump  from  s  =  0  to s  =  s, .  The  subscript  s  denotes  evaluation 
on  the  nonzero  side  of  the  shock.  The  above  equations  allow 
evaluation  of  time  derivatives  of  r„  q„  and  s,  as  algebraic 
functions  of  the  similarity  variables  and  time.  From  ( 10) 
and  ( 14), 

dV,  VA,  (,  ,  1 


From  (11)  and  (15),  we  have 
ds,  _  (  ds  \  dr], 
dt  V  dv ),  dt 


At  this  point,  the  similarity  solution  is  conceptually  com¬ 
plete.  If  one  is  given  a  physical  overdensity  field  correspond¬ 
ing  to  a  self-similar  weak  shock  at  an  instant  of  time,  r],  is  the 
radius  of  the  zero  crossing  behind  the  shock  divided  by  the 
radius  of  the  shock;  s,  is  determined  by  the  overdensity  at  the 
shock  according  to  ( 8 ),  and  t  is  the  radius  of  the  zero  cross¬ 
ing  divided  by  c0.  The  constant  of  integration  could  be  deter¬ 
mined  by  solving  ( 13)  for  a,  substituting  the  known  values 
fori,  and  t],.  From  this  point  forward,  the  shock  jump  can  be 
integrated  forward  in  time  using  ( 16).  At  any  future  time, 
the  shock  profile  as  a  function  of  range  would  be  specified  by 
(13)  from  the  zero  crossing  r  —  c^jt  out  to  the  shock  radius, 
r,  =  Cgt/rf „  with  s  =  0  for  r>r,. 

For  comparison  of  results  predicted  here  with  data,  it 
remains  to  give  an  expression  for  the  density  relaxation  time 
behind  the  shock,  r,  =  —  s(<fc/<?f)  From  equations 
(9)-(ll), 


(ds  dvV' 

r’  \dv  dt)  , 


=  -/i.  .  (i7) 

T)  ds  ,  1  +  s/2  i 

This  gives  the  density  relaxation  time,  whereas  comparisons 
are  to  be  made  with  pressure  data.  The  conversion  factor  for 
translating  density  to  pressure  relaxation  time  is  d  log/?,/ 
d  log  p„  taken  adiabatically  from  the  equation  of  state. 

The  time  coordinate  can  be  eliminated  from  ( 14)  and 
(16)  to  yield  an  equation  in  which  the  shock  radius  is  the 
independent  variable: 

dlogs,  _  r,  ds,  /dr,  =  1  -  7/(1  +  s/2) 

d  log  r,  s,  dt/  dt  1  —  i/(  1  +  s)  , 

(18) 

The  right-hand  side  of  (18)  involves  only  similarity  vari¬ 
ables,  so  that  from  (18)  and  (13)  one  can  construct  the 
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overdensity  versus  range  relation.  The  result  can  then  be 
expressed  as  overpressure  via  the  equation  of  state. 

IV.  ASYMPTOTIC  BEHAVIOR  FOR  SMALL 
OVERDENSITY 

To  demonstrate  that  (18)  leads  to  an  approximate  pow¬ 
er  law  relation,  we  truncate  ( 13)  at  first  order  in  sand  substi¬ 
tute  into  (18)  to  find 

d-°-S-  =  -  1  + - - - +  0(5,).  (19) 

d  log  r,  2  log  |os,/2|  +  4 

Anticipating  that  the  constant  a  is  of  order  unity,  the  de¬ 
nominator  on  the  right-hand  side  of  ( 19)  is  a  large  negative 
number  for  small  s,  and  the  right-hand  side  is  of  the  form 
—  1  —  6,  where  e  is  a  small,  slowly  varying  positive  quanti¬ 
ty.  Thus  an  approximate  power  relation  between  overdensity 
and  range  is  to  be  expected.  In  the  limit  of  smalls,  (19)  leads 
to 

s,  =* const  r~ 1  |log s, | "  'n.  (20) 

Since  s,  and  r,  are  approximately  inversely  proportional,  the 
logarithm  on  the  right-hand  side  of  (20)  can  be  replaced  to 
read 

s,  =  const  r,~ '  [log(r,/r0)  ]  ~  W2,  (21 ) 

where  r0  is  an  arbitrary  constant.  This  asymptotic  behavior 
recovers  a  well-known  result.311 

V.  PLANAR  FORM:  W  WAVES 

The  planar  equivalent  of  Eq.  ( 1 1 )  yields  a  solution  that 
can  be  shown  to  satisfy  the  Euler  equations  exactly,  and 
agree  with  an  assumed  adiabatic  equation  of  state  through 
second  order  in  the  overdensity.  For  a  plane  wave,  the  1/r 
term  is  absent  from  ( 9 )  and,  as  a  result  the  right-hand  side  in 
(11)  is  zero.  Removing  a  factor  ds/drj  (i.e.,  ruling  out  the 
trivial  solution  s  =  constant),  Eq.  (11)  becomes  an  algebra¬ 
ic  expression  whose  only  solution  is 

f  \ 

_  s  I 


/0=A>[1  +  (r-cj)/pc<f).  (22) 

At  any  given  time,  the  spatial  form  of  the  density  perturba¬ 
tion  is  a  linear  ramp.  When  the  ramp  has  a  zero  crossing 
midway  between  shock  discontinuities  of  opposing  sign,  the 
result  is  an  N  wave. 

The  following  procedure  can  be  used  to  show  that  (22) : 
satisfies  the  fluid  equations  to  the  stated  accuracy.  Equation 
( 22 )  for  the  density  can  be  substituted  into  the  mass  continu¬ 
ity  equation,  allowing  a  direct  integration  for  the  mass  flux 
pv.  The  mass  flux  and  density  can  then  be  substituted  into 
the  momentum  equation,  allowing  integration  for  the  pres¬ 
sure.  The  result  is  in  agreement  with  a  perfect  gas  law  with 
error  cubic  in  pc 

Shock  parameters  for  the  planar  similarity  solution  are 
easily  integrated  from  (22)  and  ( 15).  The  shock  range  r„ 
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shock  jump  s„  and  relaxation  time  r,  behind  the  planar 
shock  are 

r,  =  c0(t  +  bt'n), 

s, =bt~ l/2,  (23) 

t,  =  bt  ,n/(  1  +  bt  ~'n), 

where  b  is  an  arbitrary  constant.  In  the  planar  solution,  the 
shock  overdensity  is  predicted  to  fall  off  as  approximately 
the  inverse  square  root  of  range,  while  the  relaxation  time 
behind  the  shock  is  predicted  to  increase  approximately  with 
the  square  root  of  range. 

VI.  CONSERVATION  PROPERTIES 
A.  Mass  Integral 

Mathematical  conservation  properties  result  from  Eqs. 
(9)  and  ( 14)  for  any  solution  s(r,t).  These  properties  were 
not  maintained  nor  discussed  in  the  previous  similarity  theo¬ 
ry  for  weak  shocks.2  This  point  does  not  invalidate  that  theo¬ 
ry,  it  having  been  based  on  empirical  observation. 

In  applying  the  mathematical  solution  given  here  to 
physical  systems,  one  must  bear  in  mind  the  caveat  that  neg¬ 
ative  fluid  pressure  cannot  exist.  A  negative  density  pertur¬ 
bation  can  exist  only  to  the  extent  that  the  total  pressure 
remains  non-negative.  Density  perturbations  below  this  lim¬ 
it  in  water  would  indicate  regions  containing  cavitation  bub¬ 
bles. 

The  plane-wave  solution  (22)  and  (23)  possesses  a 
mass  integral: 

M  =  f  s(r,t)  dr  =  ( 24) 

Mass  fluctuations  ahead  of  and  behind  the  zero  crossing  at 
r  =  Cot  are  separately  conserved.  This  follows  straightfor¬ 
wardly  from  Eqs.  ( 22 )  and  ( 23 ) .  Spherical  wave  solutions  of 
Eq.  (9),  however,  do  not  explicitly  conserve  mass  for  two 
reasons.  The  first  is  the  omission  of  high-order  terms  in  ex¬ 
tracting  Eq.  (2)  from  the  Euler  equations.  The  second  and 
more  basic  reason  is  the  consistency  condition  ( 7 )  on  outgo¬ 
ing  spherical  wave  solutions  of  the  linear  wave  equation  dis¬ 
cussed  earlier.  The  spherical  analog  of  (24)  is 

M&,  =4wf  r*s(r,t)  dr,  (25) 

where  a  continuous  profile  is  assumed  between  shocks  at  r, 
and  r2.  Outside  (r,,r2),  s  is  taken  to  be  zero.  Continuous 
nonlinear  pulses  are  tactily  included  in  this  discussion  by 
allowing  that  the  magnitude  of  the  shock  jumps  at  r,  and/or 
r2  may  be  zero. 

Let  us  consider  a  solution  s(r,t)  which  may  or  may  not 
be  self-similar  and  calculate  from  (9)  and  ( 14)  the  time 
derivative  of  : 

M**  =*vrr\\r,  -!  **£  r1—  dr 

=  4ffc0J^  r+{s)  dr  =  4 irc0J  +  yj  dr.  ( 26) 

In  analogy  with  the  linear  consistency  condition  (7),  mass 
conservation  requires  the  last  integral  on  the  right-hand  side 
to  be  zero.  Physically,  if  a  wave  is  generated  that  does  not 
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satisfy  the  integral  condition,  it  will  respond  by  producing  a 
trailing  disturbance  that  offsets  the  imbalance.  The  next  sec¬ 
tion  shows  that  mass  conservation  depends  upon  impulse 
conservation,  just  as  in  the  linear  case. 

B.  Impulse  conservation 

Equation  (9 )  possesses  an  integral  of  the  motion  for  any 
solution  s{r,t),  which  is  a  generalization  of  the  linear  im- 
pluse  integral  (7).  The  result  should  be  viewed  as  an  exact 
result  of  an  approximate  equation,  so  that,  in  fact,  it  gives  a 
quantity  that  is  conserved  to  second  order  in  the  density 
perturbation.  Defining 

/=  f'  rsdr,  (27) 

Jcul 

where  r,  denotes  either  r^  or  r2  above,  Eqs.  (9)  and  (14)  give 

U'"-+Z' 

+  t)-c°ZH’+t)‘"'-0-  i28) 

In  (27)  and  (28),  it  is  not  necessary  to  specify  whether  rs  is 
greater  than  or  less  than  c0t.  The  result  is  that  integrals  ahead 
of  and  behind  r  =  c0t  are  conserved  separately.  Thus,  even 
allowing  for  the  dissipation  implicitly  present  at  a  shock 
front,  the  integral  (27)  remains  constant  for  a  spherically 
expanding  shock  to  second  order  in  the  overdensity.  Contin¬ 
uous  profiles  are  included  by  admitting  a  shock  jump  of  zero. 

With  this  result,  it  is  now  possible  to  give  the  condition 
for  mass  conservation  analogous  to  (7)  for  linear  waves.  In 
(26),  r,  and  r2  refer  to  locations  of  shock  jumps  of  opposing 
sign  with  a  continuous  profile  crossing  zero  in  between.  If 
the  integral  ( 27 )  for  the  positive  portion  of  the  profile  exact¬ 
ly  offsets  that  for  the  negative  portion,  then  the  integral  of 
the  linear  term  in  the  last  line  of  ( 26 )  is  zero  for  all  time.  The 
quadratic  term  in  (26)  must  now  be  considered.  The  limits 
of  integration  r,  and  r2  can  be  shown  to  differ  by  an  amount 
that  scales  roughly  with  the  magnitude  of  s.  This  can  be 
argued  directly  from  ( 14 ) .  [  More  convincingly  for  self-simi¬ 
lar  solutions,  the  discussion  of  Fig.  I  below  reveals  that 
r,  —  cQt^rs  (2. 5s,) 083  over  much  of  the  domain  of  interest.] 
Then,  the  integral  of  the  quadratic  term  in  the  last  line  of 
(26)  scales  roughly  as  the  cube  of  the  amplitude  of  s.  Thus, 
when 

rsdr  =  0,  (29) 

mass  is  conserved  to  within  the  accuracy  of  the  reduced  wave 
equation,  (9)  i.e.,  to  second  order  in  the  overdensity.  Eq. 
(29)  is  the  proper  generalization  of  (7).  In  fact,  they  both 
can  be  stated  in  the  form  (29),  provided  that  the  limits  of 
integration  for  linear  wave  propagation  be  interpreted  not  as 
shock  locations,  but  as  limits  outside  which  the  density  per¬ 
turbation  is  zero. 

VII.  RESULTS 

Figure  1  gives  the  similarity  solution  ( 13)  for  three  val¬ 
ues  of  the  constant  of  integration  a.  The  particular  values 
0.5,  0,7,  and  0.9  were  chosen  to  fit  relaxation  time  data  as 
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FIG.  1 .  Scaled  density  perturbation  s  as  a  function  of  t j  =  cj  fr  from  Eqs. 
(8)  and  (13).  Curves  are  labeled  with  the  values  of  a  used.  Log-tog  format 
reveals  an  approximate  power  law  relation  between  s  and  ( 1  —  7). 


discussed  below,  and  to  illustrate  the  sensitivity  of  experi¬ 
mental  measurables  to  profile  changes.  For  s  values  between 
roughly  10-4  and  10  '2,  an  approximate  power  law  relation 
exists  withr~0.4  (1  —  rj)11. 

The  overpressure  versus  range  is  integrated  numerically 
from  (18)  with  ij  given  in  terms  of  s  in  ( 13).  Two  free  pa¬ 
rameters  must  be  specified:  the  pressure  level  at  the  initial 
range  point  and  a  in  (13).  In  Fig.  2,  the  result  for  the 
overpressure  versus  range  is  given,  with  the  initial  pressure 
set  to  a  value  that  gives  reasonable  agreement  with  Eq.  ( 1 ). 
Pressures  are  calculated  from  an  equation  of  state  truncated 
after  second  order  in  j: 


P 1 

p<d 


-7+<*-«) 


(30) 


where/?  is  taken  to  be  3.5  for  water.  Shown  also  in  Fig.  2  are 
results  from  Rogers’  theory.2  Our  theory  gives  higher  pres¬ 
sures  at  close  ranges  than  Eq.  ( 1 )  This  is  reasonable  since 
overpressures  should  fall  off  faster  than  ( 1 )  at  close  ranges. 

Figure  3  gives  pressure-time  series  at  a  fixed  range  from 
Eq.  (13)  for  c0  =  1530  m/s.  These  curves  do  reflect  an  ap¬ 
proximate  exponential  decay  after  the  shock  discontinuity, 
in  agreement  with  observations.  It  is  apparent  that  the  ener¬ 
gy  and  momentum  in  the  shock  increase  as  a  a  decreases.  In 
an  actual  explosion,  the  nature  of  the  overdensity  profile 
produced  would  depend  on  the  details  of  the  energy  release 
process  and  would  likely  involve  the  Mach  number  of  the 
shock  front  while  the  shock  is  still  strong. 

In  Fig.  4,  pressure  relaxation  times  versus  range  are  giv¬ 
en  from  the  following  expression  derived  from  Eqs.  (17) and 
(30): 


t.  =  r,  JJSiL  =  r, [1  -  ( 1  -  0  ~')s)  +  0(x2). 
log  p, 

(31) 
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FIG.  2.  Shock  pressure  versus  range  from  the  empirical  power  law  ( 1 ); 
from  numerical  integration  of  (13),  (18).  and  (30)  (labels  given  in  the 
legend  are  the  a  values  used):  and  from  Rogers’  theory1:  (a)  log-log  plot; 
(b)  semilog  plot  to  magnify  differences. 


Higher-order  terms  are  estimated  to  change  this  result  by 
less  than  4%,  which  is  much  less  than  experimental  uncer¬ 
tainty.  Agreement  with  data  is  better  than  Rogers’  theory2 
predicted  at  close  ranges.  This  appears  paradoxical,  since 
Rogers’  result  for  peak  pressure  is  closer  to  the  data  in  this 
region.  Both  theories  do,  however,  tend  to  flatten  out  and 
underprcdict  the  data  at  higher  ranges.  One  likely  explana¬ 
tion  for  departure  from  data  is  that  a  real  explosion  in  gen¬ 
eral  would  not  produce  an  initially  seif-similar  profile,  but 
such  a  profile  could  evolve  over  a  finite  propagation  path. 
The  data  shown  in  Fig.  4  are  gathered  from  three  sources,  as 
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FIG.  3.  Pressure-time  series  at  fixed  range  as  computed  from  (13)  and  ( 30 ) 
for  the  specified  values  of  a,  with  p m  taken  from  the  result  of  Fig.  2. 


indicated.  It  is  unfortunate  that  no  one  data  set  covers  more 
than  a  factor  of  approximately  25  in  scaled  range.  The  simi¬ 
larity  solution  for  a  =  0.7  fits  the  Cole  data  set  reasonably, 
but  the  two  Arons  data  sets  are  not  described  particularly 
well  by  any  of  the  theoretical  models,  and  in  fact,  seem  to 
disagree  with  each  other.  Thus  the  behavior  of  the  relaxation 
time  at  long  ranges  remains  only  partially  understood.12 
This  could  indicate  mechanisms  not  included  in  the  theory. 
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FIG.  4.  Relaxation  times  r,:  experimental  values  plotted  as  individual 
points,  with  theoretical  curves  taken  from  the  empirical  power  law  (1); 
from  Eqs.  ( 17)  and  (32)  of  the  current  theory,  and  from  Rogers'  theory1 
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One  possibility  is  depth  variation  of  ambient  pressure  (or 
other  parameters)  neglected  here.  Effects  of  environmental 
variation  would  emerge  supposedly  at  long  ranges;  this  is,  in 
fact,  where  theory  and  data  begin  to  disagree  concerning  the 
behavior  of  the  relaxation  time. 

Vill.  SUMMARY 

The  nonlinear  progressive  wave  equation  (NPE)  is  an 
approximation  to  the  Euler  equations  for  nonlinear 
compressional  waves  in  an  inviscid  fluid  and,  in  fact,  is  the 
time  domain  equivalent  to  the  parabolic  equation  (PE) 
when  nonlinearity  is  omitted.  In  spherical  symmetry,  the 
NPE  possesses  a  similarity  solution  that  accounts  for  the 
following  basic  features  of  spherical  shock  waves  generated 
in  underwater  explosions:  ( 1 )  an  approximate  power  law 
relation  between  shock  overpressure  and  range;  (2)  shock 
relaxation  time  increasing  roughly  as  a  power  law  for  ap¬ 
proximately  1  decade  in  range  near  the  source;  and  ( 3 )  pres¬ 
sure-time  profiles  resembling  exponential  falloff  from  the 
shock  peak.  As  in  a  previous  semiempirical  similarity  theo¬ 
ry,2  two  free  constants  determine  the  solution.  In  the  present 
theory,  they  are  the  overpressure  at  an  initial  range  point, 
and  the  constant  of  integration  a  in  Eq.  (13).  For  plane 
waves,  the  theory  predicts  N  wave  shock  solutions  with  peak 
overpressure  varying  as  the  inverse  square  root  of  the 
elapsed  time  and  relaxation  time  varying  directly  as  the 
square  root  of  the  time.  Both  of  these  lead  to  approximate 
power  law  variation  with  range.  Only  one  free  parameter  is 
available  for  plane  waves.-  the  shock  strength  at  a  given 
range. 

For  discrete  nonlinear  pulses  with  or  without  shocks, 
the  present  theory  predicts  an  integral  of  the  motion  given  in 
Eq.  (27)  for  solutions  s(r,f).  In  Eq.  (29),  the  theory  also 
demonstrates  that  this  integral  is  a  generalization  of  the  lin¬ 
ear  impulse  integral  (7). 

Discrepancies  between  the  current  theory  and  data  are 
likely  due  to  four  assumptions  made  here  to  simplify  the 
solution.  First  is  the  reduction  of  the  fluid  equations  to  a 
first-order  nonlinear  wave  equation;  second  is  the  neglect  of 
cubic  and  higher-order  terms  in  the  density  perturbation; 
third  is  the  assumption  of  self-similarity;  and  fourth  is  the 
assumption  of  a  homogeneous  medium  ahead  of  the  shock. 
One  is  assured  that  in  the  early  phases  of  shock  wave  devel¬ 
opment  (e.g.,  in  explosive  energy  release  or  passage  of  a  su¬ 
personic  projectile)  only  the  fourth  is  reasonable.  The  degree 
to  which  the  current  theory  describes  experimental  results  in 
later  stages  is  then  a  statement  about  the  dying  out  of  tran¬ 
sients  and  the  progression  toward  a  final  state  described  by 
self-similarity. 
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